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Symbiotic nitrogen fixation (SNF) in root nodules of grain legumes such as chickpea is a 
highly complex process that drastically affects the gene expression patterns of both the 
prokaryotic as well as eukaryotic interacting cells. A successfully established symbiotic 
relationship requires mutual signaling mechanisms and a continuous adaptation of the 
metabolism of the involved cells to varying environmental conditions. Although some of 
these processes are well understood today many of the molecular mechanisms underlying 
SNF, especially in chickpea, remain unclear. Here, we reannotated our previously published 
transcriptome data generated by deepSuperSAGE (Serial Analysis of Gene Expression) to 
the recently published draft genome of chickpea to assess the root- and nodule-specific 
transcriptomes of the eukaryotic host cells. The identified gene expression patterns 
comprise up to 71 significantly differentially expressed genes and the expression of 
twenty of these was validated by quantitative real-time PCR with the tissues from five 
independent biological replicates. Many of the differentially expressed transcripts were 
found to encode proteins implicated in sugar metabolism, antioxidant defense as well 
as biotic and abiotic stress responses of the host cells, and some of them were already 
known to contribute to SNF in other legumes. The differentially expressed genes identified 
in this study represent candidates that can be used for further characterization of the 
complex molecular mechanisms underlying SNF in chickpea. 
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INTRODUCTION 

Nitrogen, besides phosphorus and potassium, is one of the pri- 
mary macronutrients for plants, and consequently a limiting 
factor for the growth of crops (Crawford, 1995). Although di- 
nitrogen (N2) represents the most abundant gas in the atmo- 
sphere, only certain prokaryotes are able to reduce it to an organic 
form (ammonia), and can thereby make it available for the growth 
of higher plants. Among these prokaryotes, rhizobia, a class 
of nitrogen-fixing bacteria, establishes symbiotic partnerships 
within the roots of some higher plants, which in turn supply them 
with energy and protect the very sensitive N2-fixation machin- 
ery from deleterious oxygen. Nutritionally and ecologically seen, 
symbiotic nitrogen fixation (SNF) identifies chickpea (and other 
legumes in general) as an important crop, and consequently a 
promising target for research (Saxena et al., 2012; Varshney et al., 
2012). Grain legumes such as chickpea form specialized organs 
in response to rhizobia invasion, the root nodules. Nodules are 
highly complex structures that protect the oxygen-sensitive N2- 
fixation machinery from its own byproducts: reactive oxygen 
and reactive nitrogen species (ROS and RNS, respectively). These 
are derived from the high rates of respiration, and the leak of 



electrons from redox proteins to O2. As a consequence, nodules 
are particularly rich not only in quantity, but also diversity of 
antioxidant defense mechanisms (reviewed in Matamoros et al., 
2003), and especially the nodule mitochondria as the principal 
source of ROS exhibit a comprehensive antioxidant repertoire 
(Iturbe-Ormaetxe et al, 2001). 

Chickpea root nodules are indeterminate in structure, since 
they maintain a persistent meristematic tissue that produces new 
cells for growth and renewal of the nodule (Lee and Copeland, 
1994). The symbiotic interaction between legumes and rhizobia 
relies on mutual signal recognition by both partners. Rhizobial 
signaling molecules called nodulation factors (Nods) are first per- 
ceived by cells in the host root epidermis, and induce the expres- 
sion of early nodulation genes (eNods) in these cells (reviewed in 
Oldroyd and Downie, 2008). Bacterial invasion of the host cells 
can occur either through root hair curls or cracks in the epi- 
dermis (Gage, 2004). The latter facilitates bacterial invasion of 
cortical cells, and does not necessarily involve Nod signaling. In 
general, the formation of indeterminate nodules is accomplished 
by root hair invasion starting with the adhesion of bacteria to root 
hairs. Subsequently, the root hairs curl, and the bacteria invade 
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the plant by a newly formed infection thread. At the same time, 
a nodule primordium is shaped by dividing cortical cells, and 
once the infection thread reaches the primordium, the bacteria 
are released into the cytoplasm of the host cells and surrounded 
by a plant-derived peribacteroid membrane, henceforth termed 
PBM (Mylona et al, 1995). PBM biogenesis and metabolism 
is governed by differential gene expression patterns of both the 
eukaryotic host legume and the prokaryotic rhizobia, which syn- 
ergistically induce the synthesis of nodulins, bacteroidines, fatty 
acids, polysaccharides, and other components. The mature PBM 
provides selectivity for metabolite and ion transport, and facili- 
tates signaling between both the prokaryotic bacteria and eukary- 
otic plant cells (Krylova et al, 2007). Although these processes are 
well understood many of the molecular mechanisms underlying 
SNF in chickpea root nodules remain unclear. 

Presently, next-generation sequencing (NGS)-coupled, 
genome-wide transcriptome profiling techniques represent the 
principal tools to interrogate the molecular mechanisms of gene 
expression in organisms across all taxa and in a wide variety of 
contexts. Especially whole transcriptome shotgun sequencing 
(RNA-Seq) has recently emerged as a potent technique for 
transcriptome studies (Libault et al., 2010; Hayashi et al., 2012; 
Reid et al, 2012; Barros De Carvalho et al., 2013). However, 
tag-based approaches such as SuperSAGE have several advantages 
compared to techniques as e.g., RNA-Seq. Transcript abundances 
are determined more accurately because of the uniform tag 
length, which impedes an introduction of biases during PCR, and 
the formation of di-tags, which allows for discrimination of PCR- 
derived tags. Additionally, the fact that only one tag is generated 
out of each transcript enables an unequivocal quantification of 
reads from a given mRNA species, and naturally results in an 
increased coverage, which facilitates the study of comprehensive 
transcriptomes as e.g. in plants (Asmann et al., 2009). In line with 
this, SuperSAGE was applied for the first time to simultaneously 
assess the differentially expressed genes from the two interacting 
organisms in Magnaporthe grisea (blast)-infected rice leaves 
(Matsumura et al., 2003). RNA-Seq, on the other hand, provides 
additional qualitative information, such as isoform expression. 
The substantially increased sequencing coverage of tag-based 
approaches must therefore carefully be weighed against the gain 
of information using whole transcriptome shotgun sequencing. 

To date, the adaption of SuperSAGE to NGS, termed deepSu- 
perSAGE, has been used to assess a broad spectrum of transcrip- 
tomes in many species (Sharbel et al., 2010; Zawada et al., 2011; 
Lenz et al., 2013) including our own works on chickpea (Molina 
et al., 2008, 2011). However, in our previous work on drought- 
and salinity-stressed chickpea roots, deepSuperSAGE transcrip- 
tion profiling was hampered by the lack of a genomic sequence of 
chickpea that prevented a faithful functional annotation of many 
SuperTags. The newly reported drafts of the chickpea genome 
(Jain et al., 2013; Varshney et al, 2013) finally changed this sit- 
uation and now enabled us to assign the majority of expressed 
SuperTags in roots and nodules to a genomic locus and thereby to 
a potential function in the context of nodulation. The approx- 
imately 50,000 duplicate and homopolymer-filtered SuperTags 
from our previous study in fact represent nearly 1800 genes of 
which at least 800 are commonly expressed in both tissues. Up 



to 682 are more abundant in nodule tissue, and 71 genes dis- 
play a highly significant differential expression (p-value < 0.01). 
The underlying data integrity was previously confirmed via 
microarray hybridization of approximately 3000 UniTags with 
diverse regulation levels. Of these, 660 could be reliably measured 
via hybridization, and the comparison between both platforms 
resulted in a shared tendency toward up- or downregulation of 
these transcripts of 79% (Molina et al., 2011). The identified set 
of differentially expressed genes consequently reflects necessary 
adaptions of the host cell transcriptome with respect to SNF in 
chickpea nodules. 

In the past, many aspects of nodule development in legumes 
have been thoroughly characterized (see Ferguson et al., 2010; 
Desbrosses and Stougaard, 2011; Hayashi et al., 2013). However, 
less emphasis was paid to the transcriptomes of both nodules and 
roots, especially in chickpea, since annotation had to be based 
on the genome sequences of other legumes (e.g., Medicago truti- 
catula). Now that the draft genome sequence is available, we 
reanalyzed the transcriptomes of unstressed chickpea nodule-free 
roots and mature root nodules from our previous study, and 
confirmed the newly identified expression patterns via quanti- 
tative real-time PCR (qRT-PCR) using five individual biological 
replicates. 

MATERIALS AND METHODS 
PLANT MATERIAL 

Beja 1 (INRAT 93-1) is a salt-tolerant chickpea (C. arietinum) 
variety (L'Taief et al., 2007) that was released in 2003 from the 
National Institute for Agricultural Research of Tunisia (INRAT) 
and the International Center for Agricultural Research in the Dry 
Areas (ICARDA). Rhizobial inoculations and growth conditions 
for chickpea plants are described in the work of Molina et al. 
(2011). Briefly, surface-sterilized Beja 1 seeds were germinated on 
0.9% agar for five days in a dark chamber at room temperature, 
and seedlings with a minimum root length of 5 cm were sub- 
sequently inoculated with Mesorhizobium ciceri strain UPMCa7 
(Romdhane et al., 2007). The inoculated seedlings were hydroare- 
oponically grown in a temperature-controlled glasshouse with a 
day/night temperature regime of 28/20°C, respectively, and a 16 h 
photoperiod for 40 days. After this period, the six-week-old chick- 
pea plants were harvested from the hydroponic cultures. Nodule 
tissue was carefully separated from the remaining root system, 
and subsequently both tissues were immediately stored in liquid 
nitrogen. Plant breeding was carried out in the greenhouse facil- 
ities of the "Soil and Symbiosis Research Unit" of the National 
Institute for Agricultural Research (INRA) in Montpellier, France. 

TOTAL RNA ISOLATION, LIBRARY PREPARATION AND 454 SEQUENCING 

Dissected mature nodules were used in their entirety for charac- 
terization of the nodule-specific transcriptome, while transcrip- 
tion profiling of the roots was performed with all the remaining 
root material available. Total RNA isolation, and subsequent 
construction and sequencing of SuperSAGE libraries from these 
tissues were performed as previously described using RNA from a 
pool of 15 plants (Molina et al, 2011). In brief, total RNA was 
isolated as described by Pawlowski et al. (1994) with a modi- 
fied precipitation of the RNA in 3M LiCl at 4°C overnight. Then 
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the polyadenylated fraction of the total RNA was purified using 
the Oligotex mRNA Mini Kit (QIAGEN, Hilden, Germany), and 
subsequently used for construction of SuperSAGE libraries as 
detailed by Matsumura et al. (2003) with a modified sequenc- 
ing procedure. Instead of di-tag concatenation and subsequent 
cloning for Sanger sequencing, the PCR amplified di-tags were 
directly sequenced on the GS20 platform (454 Life Sciences, 
Branford, USA). 

BIOINFORMATICAL PROCESSING OF deepSuperSAGE SEQUENCING 
DATA 

Sequencing data was analyzed with the GenXPro SuperSAGE 
data processing pipeline. First, distinct libraries were sorted out 
from the bulk of sequences according to their respective bar- 
codes. Then, PCR-derived and all low-complexity reads con- 
taining 12 or more consecutive adenine bases were eliminated. 
These filtered tags were subsequently mapped against the recently 
published draft of the chickpea genome (Varshney et al., 2013) 
using the short read mapper Novoalign v2.07.13 (Novocraft 
Technologies) with default parameter settings. Tags mapping to 
more than one locus were excluded from further analysis. Finally, 
feature annotation for the mapped loci was performed with 
the mRNA sequences from the "Official Gene Set" (OGSvl.O; 
Varshney et al, 2013), and reads were counted using the Python 
package HTSeq v0.5.4p2 (EMBL Heidelberg, https://pypi. python. 
org/pypi/HTSeq). The numbers of unambiguously annotated 
SuperTags were normalized to 10,000 sequenced tags in total (tags 
per ten thousand; TPT) for each library to warrant comparabil- 
ity between the libraries. Statistical significance was assessed by x 2 
tests (Man et al., 2000), and fold changes were determined by pair- 
wise comparison of the normalized tag numbers. TPT counts of 
zero were adjusted to 0.05 to allow for calculation of fold changes 
even if a given tag was only present in one of the libraries (see 
Table SI). 

Functional annotation of the expressed genes was performed 
with the MapMan software (version 3.5.1) developed by Thimm 
et al. (2004). First, a reference mapping file comprising a draft 
metabolic network of chickpea was generated by classification 
of all protein-coding sequences from the draft genome into 
MapMan functional categories via the Mercator tool (Lohse 
et al, 2014). Approximately 60% of these sequences could be 
assigned to one of the 34 functional classes (Figure SI, Table S2). 
Subsequently, the generated reference file was used for mapping of 
the differentially expressed genes onto different metabolic path- 
ways and to assign these genes to several large enzyme families. 

CONFIRMATION OF THE GENOME-BASED REANALYSIS BY 
QUANTITATIVE REAL-TIME PCR 

Quantification of mRNA by real-time PCR was performed on 
the StepOne Real-Time PCR System (Applied Biosystems) using 
independent biological replicates that were bred in the same way 
as described above. Root and nodule tissues from 6 freshly grown 
plants at the age of six weeks were dissected and used for total 
RNA isolation with the InviTrap Spin Plant RNA Mini Kit (Stratec 
Biomedical) following the recommendations of the manufacturer 
for use of lysis solution DCT While nodules were simply stripped 
off the snap-frozen material and afterwards used entirely for total 



RNA isolation, dissection of root tissue was performed with a 
sterile razor to obtain absolutely nodule-free root tissue for char- 
acterization of the root-specific transcriptomes. Remaining DNA 
fragments in the isolated total RNA were digested by DNase I 
(Baseline-ZERO, Epicentre) as recommended by the manufac- 
turer. Subsequent to quantification of the total RNA (Qubit, 
Life Technologies), all isolates were quality-controlled on the 
Bioanalyzer 2100 (Agilent Technologies). Isolates with an RNA 
Integrity Number (RIN) of 7 or higher were reverse-transcribed 
with Superscript III Reverse Transcriptase (Invitrogen) following 
the manufacturer's instructions for first-strand cDNA synthesis. 
Reverse-transcribed cDNA corresponding to 20 ng total RNA was 
then added to each amplification reaction on the StepOne Real- 
Time PCR System. All amplification reactions were carried out 
in 12 ul volume with the 5x HOT MOLPol EvaGreen qPCR Mix 
(ROX) from Molegene, complemented by the respective forward 
and reverse primers (Table 1) in a final concentration of 250 nM 
each. Initial denaturation was performed at 95° C for 15 min, fol- 
lowed by 40 cycles of 15 s at 95° C, 20 s at 65° C and 30 s at 72° C. 
A final elongation step at 72° C for 5 min allowed the polymerase 
to complete all unfinished strands. Subsequently, a melting curve 
analysis was performed to verify exclusive amplification of the 
expected products. Additionally, the threshold cycles (C f values) 
of the negative (no template) controls from all employed assays 
were ensured to be higher than 35. 

The relative transcript abundances between root and nod- 
ule tissue from the different biological replicates were calculated 
according to the AAQ method using the geometric mean of 
three previously determined reference genes. All target and ref- 
erence mRNAs were quantified in duplicates for all biological 
replicates, respectively. The arithmetic mean of each duplicate 
was then used to calculate the AAC f values between the tis- 
sues. Thirteen candidate reference genes were screened for their 
target stability using a pool of reverse-transcribed cDNAs from 
the five biological replicates that passed quality control (Figure 
S2A). The determined expression ratios were analyzed with 
geNorm (Vandesompele et al., 2002), and six of the best per- 
forming candidates (geNorm M < 0.5) were analyzed in more 
detail using three individual biological replicates (Figure S2B). 
The optimal number of target reference genes was determined 
to vary around two to three. To ensure optimal comparability, 
the 3 best-performing candidates from the individual test run 
(Ca_04993, Ca_09743, and Ca_03068) were subsequently used 
for normalization of the gene expression ratios from root and 
nodule tissue. 

RESULTS AND DISCUSSION 
GENOME BASED REANALYSIS 

With a duplicate- and homopolymer-filtered read number of 
25,160 (roots) and 26,380 (nodules), both libraries contain a 
similar number of reads (±5%). Genomic mapping resulted in 
17,909 (71.1%) and 20,508 (77.7%) unambiguously assigned tags 
(root and nodule, respectively). The remaining reads were either 
mapped to more than one locus or could not be mapped to 
the mRNA-encoding sequences of the chickpea genome at all 
(Figure 1). Only uniquely assigned reads were taken into consid- 
eration for further analysis. 
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Table 1 | List of targeted mRNAs along with the respective primer and probe sequences used for qRT-PCR quantification. 



Gene ID* 


Alias* 


Accession number 


Primer 


Sequence 


Amplicon size 


Ca_04993+ 


Tubulin alpha-7 chain 


XM. 


.004501016.1 


Forward 


GTG GTG ATCTTG CCAAG GTTCAG 


222 












Reverse 


GACTCAGCACCAACCTCTTCATAATC 




Ca. 


.09743+ 


Heat shock protein Hsp90 


XM. 


004491473.1 


Forward 


TGTTG AAG CTTG GACTG AG CATTG 


111 








XM. 


004491474.1 


Reverse 


tp/^ Ar* ot on — ro o a Ton — r o ot a o o 

TCGACCTCTTCCATCTTGCTACC 




Ca_03068+ 


S-adenosylmethionine 


XM. 


.003609813.1 


Forward 


CCTCACTATCGTGAAGAACAGCTTTG 


166 






synthase 






Reverse 


o o o a ~r — l — r a a^apo on — r o a o o a ot. — r o 

CCCATTTAAGAGGCTTCACCACTTC 




Ca. 


.06305 


Monothiol glutaredoxin-S17 


XM. 


.004504992.1 


Forward 


TGCTCCAAG ATG CG G CTTTAG 


187 












Reverse 


CCATAACAATATCGCAACCGCCTATC 




Ca. 


.13049 


Integrator complex subunit 


XM. 


.004512818.1 


Forward 


AAATTG C AG CTGATTTG G CTTC C 


302 












Reverse 


AAG CTTTGTAAG G GTCCTCTGTATG 




Ca. 


.22734 


Putative uncharacterized 


XM. 


.004512804.1 


Forward 


CAATGAAGCGTTCG G GTTTGTG 


114 






protein 






Reverse 


CTCCGACCGCCACAACATATC 




Ca. 


.10340 


Multidrug resistance protein 


XM. 


.004503208.1 


Forward 


AGAGTCAGGGCATGACACTCATC 


148 








AB024992.1 


Reverse 


CCGTGCCATGGGATGCTTAG 




Ca. 


.04229 


Neutral alpha-glucosidase 


XM. 


_UU4bUzyzo.1 


Forward 


GGCACCTACTTCTGGTGGAAATG 


169 






AB 






Reverse 


AAG CTG GTGAATGTG CCCTTTG 




Ca. 


.07680 


Putative uncharacterized 


XM. 


_UU44y4yb/. 1 


Forward 


C AG G AAAC AG CC G AAATCTAG G ATG 


132 






protein 






Reverse 


AC AAG CTTCTG G C C AACTATTG C 




Ca. 


.06862 


Putative uncharacterized 


XM. 


.004508211.1 


Forward 


G CTGTCCATGAGAAAG GAG ATGTG 


142 






protein 






Reverse 


a o ot o ot on — ro o a a a ot 'o on — l — ro 

AG CTG CTCTTGG AAACTG CTTTG 




Ca. 


_1 1013 


Putative uncharacterized 


A I VI. 


nn/i a no a o a 1 
_UU44yo4y4. I 


Forward 


ATG GTG C AC ATG G AG AATTC ATG G 


242 






protein 






Reverse 


O O A AOATAOOA AO O O OTO O AT A f* 

GLAALAI AGGAAGLLL 1 GCAI AG 




ua_ 


_ iUo I z 


Coronatine-insensitive 1 


AlVI. 


_UU4oUo I /o. I 


Forward 


a oootatootoo at ot o o a toto 

AGGGTATGGTGCATCTCCATCTG 


181 












Reverse 


AATCTGATCTTTGGCCAGCAAGAG 




La. 


1 COO A 

_ lobo4 


HMG l/Y like protein 


XM. 


.00451 2oo9.1 


Forward 


AACAACACCTGCTAGTGCTCAAC 


110 












Reverse 


AAATG AG G C CTAAG C ACTG C AAG 




Ca. 


_05800 


6-phosphogluconate 


XM. 


.004503591.1 


Forward 


CTTGTTCAGGCTCAGAGGGATTTG 


126 






dehydrogenase 






Reverse 


AATTAAGAGCAGCAACACCAGTACC 




Ca. 


o o r\ o o 

.22023 


Monosaccharid transport 


XM. 


.004504903.1 


Forward 


CATGTTGCCTGAGACTAAGGGAATAC 


133 






protein 






Reverse 


ta aoaootoo r*i — ro o o o a toto 

TAAC AG CTCCCTTGCCCATCTC 




La. 


.00007 


Squamous cell carcinoma 


XM. 


.004485354.1 


Forward 


~T — 1 — r~0 O A O O A TO A O O AO On — r~0 n — l~0 

TTTGGACGATGAGCACCTTGTTG 


225 






antigen 


XM. 


.004485355.1 


Reverse 


A ATOOTOTOAOn — TOOTOOOT — 1 — ro 

AAI GL 1 L 1 LAL 1 1 Lb 1 GGL MIL 










XM. 


.004485356.1 














XM. 


.004485357.1 








Ca. 


.05370 


Prefoldin subunit 


XM. 


.004497178.1 


Forward 


CTCAACACGTTCTCGTCGATGTC 


134 








XM. 


.004497179.1 


Reverse 


TTTG G G ATG C C AC CTC AAC AAG 










XM. 


.004497180.1 








Ca. 


.15777 


Serine/threonine protein 


XM. 


.004506314.1 


Forward 


TGGACCCTGAATACTATAGGCTACAAC 


332 






i ■ i ■ i j. f - n a 

kinase-hke protein CCR4 






Reverse 


ACGCCGCAAGAGCTGTTTC 




Ca. 


.12354 


ADP-ribosylation factor 


XM. 


.004509671.1 


Forward 


TTCCATCTCCAGTGCCGATCTC 


200 






GTPase-activating protein 






Reverse 


CAG AG AATTCG GTCTTG AAG ATCTGTC 




Ca. 


.15466 


Nodulin 6 


XM. 


.0044979371 


Forward 


ACTGATGCCTATGCATTTCCTGAAC 


124 












Reverse 


ACTTCCACAGCCTCCGGAAC 




Ca. 


.12714 


Putative uncharacterized 


XM. 


004502350.1 


Forward 


G G CC AATC CTG AG AAG AG AATC AC 


229 






protein 


XM. 


.004502351.1 


Reverse 


CCATGCCTCCTCCAACAAATTGTC 




Ca. 


.13139 


Putative uncharacterized 


XM. 


.004498271.1 


Forward 


TGG CTG AACAAACTCATTTG G GAAG 


219 






protein 


XM. 


.004498272.1 


Reverse 


CCTGCAACCTTGATATCTCCAGGAAC 




Ca. 


.03442 


Glutathione S-transferase 


XM. 


.004495920.1 


Forward 


GGAAGAGAATGAAGCCAAGTTGAACAC 


217 












Reverse 


TAGACCAAG CTG GTCTTG CAGTG 




Ca. 


.16084 


Leghemoglobin 


XM. 


004490852.1 


Forward 


GAGATG CTACATTG GGTG CTGTTC 


159 












Reverse 


GCCAATCCATCATAGGCGAGTTC 





*Gene ID and alias according to OGSvl.O (Varshney et al., 2013); f NCBI reference sequences for all transcript variants targeted by the respective primer pair; 
+ reference gene used for normalization. 
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In our study in 2011, the 51,545 sequenced Beja 1 SuperTags 
from root and nodule tissue of unstressed chickpea plants were 
grouped into 11,525 UniTags (SuperTags of common origin), 
while the genome-based reanalysis of these SuperTags resulted 
in the identification of approximately 1780 expressed gene loci. 
These represent the transcriptionally most active genes in chick- 
pea roots and nodules, and consequently the corresponding 
mRNAs are the most abundant transcripts. The 140 UniTags pre- 
viously found to be more than 8-fold prevalent in nodules were 
reduced to 61 significantly differentially expressed loci in the 
genome of chickpea (Table 2). On the other hand, the number of 
more than 20-fold differentially expressed UniTags between both 
tissues increased from four transcripts to 64 that actually exhibit 
more than a 128-fold differential expression. Consequently, many 
of the slightly upregulated UniTags in nodules are expressed from 
the same genomic locus, and based on the genome sequence these 
can be combined to provide more meaningful data. 

NODULE-SPECIFIC GENE EXPRESSION OF CHICKPEA CULTIVAR Beja 1 

The identified expression profiles of Beja 1 root and nodule tis- 
sues share around 800 (~45%) expressed genes (Figure 2), while 
682 (almost 40%) of the genes are heavily upregulated in nodule 



tissue (vs. 297 distinctly expressed genes in roots). Although the 
captured number of expressed genes (~1780) is relatively low 
compared to more recent profiling studies, the large difference 
in distinctly expressed genes indicates important variations in 
expression of the most abundant transcripts in the context of 
nodulation. The relatively high number of expressed genes in 
nodule tissue reflects an induced expression of a plethora of genes 
that are putatively involved in establishment and maintenance of 
the symbiotic relationship. A total of 71 genes were identified 
as highly significant (a = 0.01) differentially expressed between 
both tissues (Table SI). The distribution of these genes across 
the 8 chickpea chromosomes is relatively uniform and varies 
around nine (±50%) differentially expressed genes per chromo- 
some. With respect to the varying chromosome sizes, most of 
the differentially expressed genes are located on chromosome 8, 
but no particular chromosomal region was found to be signifi- 
cantly enriched with differentially expressed loci. Genes showing 
more than a 2-fold differential expression between roots and nod- 
ules are listed in Table 2. Consistent with the high number of 
expressed genes in nodule tissue, the number of upregulated genes 
in nodules is about twice as high as in roots. As expected, most 
genes display a relatively low differential expression, and there- 
fore the respective significance levels vary accordingly. Sixty four 
significantly differentially expressed transcripts (enrichment of 
128-fold or more) could be identified as highly enriched in one of 
the tissues. Among others, these genes include BZIP transcription 
factor 2 (inferred from Phaseolus vulgaris), nodulin 6, abscisic 
acid receptor PYL4, and glutathione S-transferase (inferred from 
M. truncatula) all of which are upregulated in nodules (Table SI). 
The complete set of expressed genes is depicted in a heat map 
(Figure S3) that illustrates the relatively high expression ratio in 
nodule compared to root tissue, since numerous transcripts are 
found to be more abundant in nodules. 




FIGURE 2 | Tissue-specific gene expression in chickpea cultivar Beja 1. 

Outer numbers represent genes that are expressed in root (left) or nodule 
(right) tissue, while the number within the overlap represents commonly 
expressed genes. 
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FIGURE 1 | Genomic mapping results of the captured transcriptomes 
from chickpea root and nodule tissue. 



Table 2 | Number of differentially expressed genes in root and nodule tissue from chickpea cultivar Beja 1. 

Differentially expressed 

>2-fold >4-fold >8-fold >64-fold >128-fold >256-fold 

Nodule upregulated 108(953) 92(755) 61(692) 51(296) 51(51) 5(5) 

Nodule downregulated 56(228) 42(159) 26(135) 13(50) 13(13) 3(3) 

Sum of differentially expressed genes 164(1181) 134(914) 87(827) 64(346) 64(64) 8(8) 

Numbers on the left represent significantly differentially expressed genes (a = 0.05), while the numbers in brackets include all expressed genes regardless of any 
significance threshold. 
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VALIDATION OF THE IDENTIFIED GENE EXPRESSION PATTERNS BY 
QUANTITATIVE REAL-TIME PCR 

Differential expression of the ten most up and downregulated 
genes in chickpea nodule tissue from five biological replicates 



was validated via qRT-PCR and revealed a comprehensive biolog- 
ical variance between the different isolates (Figure 3A, Table SI). 
While all of the upregulated and seven out of the ten most down- 
regulated genes in nodule tissue are found to be accordingly 
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FIGURE 3 | Biological variance in gene expression of 21 candidate and 
three reference genes across five biological replicates (A) and 
comparison of deepSuperSAGE expression patterns from ten pooled 
plants with the individual expression ratios determined by qRT-PCR (B). 

The logarithmized (base 2} expression ratios between root and nodule tissue 
from five biological replicates are depicted in box plots for the indicated 
genes. Positive values represent nodule upregulated mRNAs and negative 



values mRNAs that are more abundant in root tissue. The gene expression 
ratios of the individual biological replicates (qPCRI-5) in comparison to the 
pooled (SSage) tissues are additionally shown in the heat map. Nodule 
upregulated transcripts are represented in red and downregulated transcripts 
in green. Clustering was performed with the MultiExperiment Viewer version 
4.9 by hierarchical clustering of all genes and samples using Euclidean 
distance calculations. 
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expressed in at least one of the replicates, some of these genes 
exhibit converse expression ratios in the other replicates. The 
expression of leghemoglobin, which is known to be expressed 
by legumes in response to colonization of the roots by rhizobia 
(Benedito et al., 2008; Libault et al., 2010), was assessed addition- 
ally to the twenty most differentially expressed genes. As expected, 
the mRNA encoding leghemoglobin is significantly more abun- 
dant in nodule tissue regardless of the biological replicate. The 
individual expression ratios, however, range from almost 25,000- 
fold down to 10-fold upregulated in the respective tissues. The 
extensive biological variance in expression of some of the twenty 
candidate genes becomes even more apparent after hierarchical 
clustering (Figure 3B). The individual replicates can be broadly 
classified into two groups that show similar expression ratios 
either for the upregulated (qPCRl and qPCR2) or alternatively for 
the downregulated transcripts (qPCR3-5) compared to the pool 
of ten plants used for deepSuperSAGE. 

UPREGULATED GENES IN CHICKPEA NODULES 

The 10 most upregulated genes in chickpea nodule tissue are 
listed in Table 3. The proteins encoded by two of the listed genes 
have a putative function, but the respective protein sequences 
could be matched with the InterPro database, which pre- 
dicts a serine/threonine-protein kinase and a drug/metabolite 
transporter. Especially the transcript encoding the putative 
drug/metabolite transporter is highly upregulated in nodules 
(> 300-fold enrichment) not only in the pooled plant tissue 
but also in the individual biological replicates, which suggests 
that the corresponding gene product is functionally important 
in the rhizobia-adapted metabolism of nodules. The putative 
drug/metabolite transporter (IPR000620) is predicted to be an 
integral membrane protein and may contribute to the selectivity 
of metabolite transport through the mature PBM. 

Another strongly upregulated transcript in nodule tissue 
encodes a glutathione S-transferase (GST) that is implicated 



in the antioxidant defense of legume root nodules. GSTs can 
directly scavenge peroxides with the help of glutathione as elec- 
tron acceptor, and furthermore detoxify endogenous compounds 
such as peroxidized lipids via conjugation of glutathione to tar- 
get substrates, which facilitates their sequestration and removal 
(reviewed in Becana et al, 2010). Not surprisingly, GSTs consti- 
tute a large gene family in legumes, whose N2 -fixation machinery 
needs protection from ROS and RNS. Nodulin 6, encoded by early 
nodulin gene MtN6 in M. truncatula, is one of the early response 
genes upon rhizobia infection in legumes. The mRNAs coding 
for Nodulin 6 accumulate in outer cortical cells that contain 
pre-infection threads, and in front of growing infection threads 
in mature root nodules (Mathis et al., 1999). Upregulation of 
Nodulin 6 consequently reflects the ongoing root hair invasion 
in mature nodule tissue. 

The transcript encoding 6-phosphogluconate dehydrogenase 
(6PGD) is more than 200-fold upregulated in the pooled nod- 
ule compared to root tissue. The encoded enzyme belongs to the 
family of oxidoreductases and is involved in the non-oxidative 
phase of the pentose phosphate pathway, where it catalyzes the 
conversion of 6-phosphogluconate to ribulose 5-phosphate. The 
upregulation of 6PGD resonates with the induced expression 
of ribulose-phosphate 3-epimerase (~170-fold upregulated in 
nodules, Table SI), which in turn uses ribulose 5-phosphate as 
a substrate to generate the ketose sugar xylulose 5-phosphate. 
A primary end product of the pentose phosphate pathway is 
NADPH, which is needed in response to oxidative stress, besides 
being necessary for fatty acid synthesis. NADPH can serve as 
co-substrate for glutathione reductases that reduce oxidized glu- 
tathione e.g. subsequent to its oxidation via GSTs (Table 3) or 
other glutathione peroxidases (up to 3-fold upregulated in nod- 
ules, Table SI). Further end products of the pathway include 
ribose-5-phosphate, used in the synthesis of nucleic acids, and 
erythrose-4-phosphate, implicated in the synthesis of aromatic 
amino acids. Against this backdrop, the induced expression of 



Table 3 | Proteins encoded by the 10 most upregulated genes from nodule compared to root tissue of chickpea cultivar Beja 1. 



ID 




TrEMBL database 


InterProScan 


Fold change 


P-value 


Chr. 


Ca_ 


.13139 


Putative uncharacterized protein 


Drug/metabolite transporter 


346 


0.0008 


4 


Ca. 


03442 


Glutathione S-transferase 


Glutathione S-transferase 


346 


0.0008 


4 


Ca. 


.12714 


Putative uncharacterized protein 


Protein kinase, catalytic domain; 
Serine/threonine-protein kinase domain 


311 


0.0015 


5 


Ca_ 


.12354 


ADP-ribosylation factor GTPase-activating 
protein AGD10 


Arf GTPase activating protein 


276 


0.0028 


7 


Ca_ 


.15466 


Nodulin 6 


Amidohydrolase 2 


276 


0.0028 


4 


Ca_ 


.15777 


Serine/threonine protein kinase-like protein 
CCR4 


Protein kinase, catalytic domain; 
Serine/threonine-protein kinase domain 


242 


0.0051 


6 


Ca_ 


.05370 


Prefoldin subunit 


Prefoldin subunit 


242 


0.0051 


4 


Ca_ 


.05800 


6-phosphogluconate dehydrogenase, 
decarboxyiating 


6-phosphogluconate dehydrogenase, 
C-terminal; NAD-binding 


207 


0.0095 


6 


Ca_ 


.22023 


Monosaccharide transport protein 


Sugar/inositol transporter 


207 


0.0095 


6 


Ca_ 


.00007 


Squamous cell carcinoma antigen recognized 
by T-cells, putative 


RNA recognition motif domain; 
RNA-processing protein, HAT helix 


207 


0.0095 


1 



Protein functions inferred from the TrEMBL and InterPro database are listed along with the respective fold changes, p-values (j 2 ) and chromosomal positions. This 
list represents an excerpt from Table S 1. 
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several members of the pentose phosphate pathway is not exclu- 
sively linked to an increased anabolism in nodule tissue, but 
rather ensures a steady supply with reducing equivalents (in form 
of NADPH). 

The heavily upregulated monosaccharide transport protein 
(MTP) represents an ortholog of the hexose transporter Mtstl 
in M. truncatula. An alignment of the mRNA sequence of Mtstl 
from M. truncatula with the genomic sequence of MTP resulted 
in 44% matching bases, although the MTP sequence still retains 
introns (Figure S5). Interestingly, expression of Mtstl was asso- 
ciated with a successfully established symbiosis of M. truncatula 
and the vesicular-arbuscular mycorrhizal fungus Glomus versi- 
forme (Harrison, 1996). In situ hybridization revealed high Mtstl 
mRNA levels in phloem fiber cells of the vascular tissue, cells of 
the root tip, and in cortical cells of the mycorrhizal root, espe- 
cially in highly invaded areas. The author therefore suggests that 
increased expression levels of Mtstl are correlated with inter- 
nal growth of the fungus and with a functioning symbiosis, 
because the affected cells exhibit an increased metabolism, which 
in turn requires an intensified energy supply. To our knowledge 
an ortholog of Mtstl in chickpea has not yet been described. 
Analogous to the function of Mtstl in the context of vesicular- 
arbuscular mycorrhizal associations, MTP potentially ensures the 
sugar supply for root cells directly involved in the symbiotic 
association of rhizobia and legumes. The strong upregulation 
of MTP in nodules might furthermore contribute to the salt- 
tolerant trait of chickpea cultivar Beja 1, since hexose sugars 
also increase the cytosolic solute concentration, and thus can 



contribute to maintain the osmotic pressure (Hasegawa et al., 
2000). One of the major bottlenecks of SNF in plants is the sensi- 
tivity of the symbiotic interaction, which renders the nodules very 
susceptible to abiotic stresses. The activity of enzymes directly 
involved in SNF, for example, was drastically decreased in salt- 
stressed nodules (Cordovilla et al., 1994). The nodule-specific 
induction of MTP along with 6PGD thus indicates that several 
of the upregulated genes in nodules are not only linked to an 
increased metabolism, but also involved in the salt-tolerant trait 
of Beja 1. This is in line with our previous finding that several 
relatively low expressed UniTags in unstressed nodules become 
highly abundant in salt-stressed nodule tissue (Molina et al., 
2011). 6PGD was additionally identified as a stress-responsive 
protein in a comparative proteomic analysis of A. thaliana roots 
that had been exposed to 150 mM NaCl (Jiang et al., 2007). 
The differential expression of this gene with the onset of salt 
stress in A. thaliana seems to be tightly linked to the accumula- 
tion of compatible osmolytes, and the present findings empha- 
size the importance of this mechanism in the context of SNF 
in legumes. 

FUNCTIONAL CLASSIFICATION OF DIFFERENTIALLY EXPRESSED 
GENES IN CHICKPEA NODULE TISSUE 

The identified expression patterns of chickpea root and nodule 
tissues display comprehensive differences. Many gene products 
involved in sugar metabolism, antioxidant defense as well as 
biotic and abiotic stress responses of the host cells show dras- 
tically altered levels after adaption to the symbiotic relationship 
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FIGURE 4 | MapMan-based classification of the expressed genes into 
large enzyme families. Each field represents the expression of a particular 
gene. Upregulated genes in nodule tissue are shown in red, genes with 



reduced expression in green, and undifferentially expressed ones in white. 
Dark gray fields indicate that none of the expressed genes could be assigned 
to the respective class. 
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(Figure S4). This is in line with a recent report of nodulation- 
relevant genes identified in soybean roots 10 days after inocula- 
tion with Bradyrhizobium japonicum (Barros De Carvalho et al., 
2013). The upregulated genes in nodules were found to be pri- 
marily involved in host cell metabolism, cell wall modifications 
and the antioxidant defense system. Interestingly, glycolysis and 
the citric acid cycle were the most active metabolic pathways in 
the context of SNF. The strong induction of several members of 
the pentose phosphate shunt in Beja 1 root nodules, however, 
reflects a shifted nitrogen balance in nodule tissue that is linked to 
an increased anabolism. Naturally, the interacting host cells have 
to ensure a steady supply of metabolites for the nitrogen-fixing 
rhizobia on the one hand, but also provide nitrogenous com- 
pounds for the rest of the plant on the other. Accordingly, the gene 
encoding asparagine synthetase is heavily upregulated in nod- 
ule tissue (~70-fold, see Table SI, Figure S4B), since asparagine 
represents one of the primary nitrogen transport products in 
plants. 

The comprehensive expression of genes encoding GSTs in 
nodule tissue is accompanied by an upregulation of several oxi- 
doreductases involved in different metabolic pathways as well 
as the two isoflavone 7-O-methyltransferase isoforms 8 and 9 
(Figure 4). A globally induced expression of genes encoding GSTs 
was previously identified in nodule tissue from M. truncatula 
by Benedito et al. (2008), and the present findings confirm the 
importance of ROS scavenging for SNF in chickpea. Interestingly, 
overexpression of isoflavone 7-O-methyltransferase 8 was shown 
to enhance disease resistance in M. sativa (He and Dixon, 2000), 
while the other isoform was identified as drought-responsive 
transcript in chickpea (Varshney et al, 2009). Against this back- 
drop, upregulation of isoform 9 could well contribute to the 



salt-tolerant trait of Beja 1, since an a priori increased expression 
of stress-responsive genes already prepares the very sensitive inter- 
action of eukaryotic and prokaryotic cells in root nodules for a 
possible onset of stress. Expression profiles of other large enzyme 
families such as nitrilases, glucosidases or the cytochrome P450 
superfamily of monooxygenases are less consistent, but several 
peroxidases implicated in the response to oxidative stress as well 
as defense responses to fungi are heavily upregulated in nod- 
ules. These peroxidases apparently add to the comprehensive 
antioxidant defense repertoire that protects the oxygen-sensitive 
N2 -fixation machinery in chickpea nodules. 

Functional annotation onto known biotic stress pathways 
reveals a reduced expression of disease resistance proteins, defense 
genes and of genes involved in mediating host cell responses to 
pathogens (Figure 5). This is consistent with the massive down- 
regulation of the multidrug resistance and disease resistance 
response proteins in mature Beja 1 nodules, and confirms our 
previous observation that a sustained inhibition of at least some 
of the defense reactions in root nodules seems to be a prerequisite 
for maintenance of the symbiosis. 

The functional analysis of genes involved in generation of 
secondary metabolites in nodules illustrates a strong upregu- 
lation of genes implicated in phenylpropanoid and flavonoid 
biosynthesis although several genes encoding proteins of the 
shikimic acid pathway are downregulated (Figure 6). Especially 
dihydroflavonols appear to be of special importance for SNF, 
which is likely linked to their function as antioxidants and radical 
scavengers (Haraguchi et al., 1996). The most upregulated gene 
involved in the synthesis of phenylpropanoids actually encodes 
a nicotinamidase that operates in the salvage pathway of NAD 
biosynthesis (Wang and Pichersky, 2007). The induced expression 




FIGURE 5 I Gene expression relating to biotic stress pathways based on the functional annotation with MapMan. Please consult Figure 4 for further 
details. 
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FIGURE 6 | Differential expression of genes involved in biosynthesis of secondary metabolites. Please consult Figure 4 for further details. 



of this gene in nodule tissue underlines the importance of NAD 
as redox carrier in the rhizobia-adapted metabolism. 

CONCLUSIONS 

The increasing throughput of next-generation sequencing tech- 
nologies and the application of 3rd generation sequencing (such 
as the SMRT sequencing system from Pacific Biosciences) are 
currently intensifying the efforts to decipher the genome of a 
steadily increasing number of legumes. As in the case of chick- 
pea, these genome sequences will allow for more detailed analyses 
of the differential gene expression in relation to SNF in legumes. 
Using the recently published chickpea genome sequence, we rean- 
alyzed the transcriptomes of unstressed root and nodule tissue 
and identified a strong upregulation of genes encoding glu- 
tathione S-transferases or with implications in phenylpropanoid 
and flavonoid biosynthesis. The expression of twenty candidate 
genes from the reanalyzed dataset was validated using five indi- 
vidual biological replicates, revealing a comprehensive biological 
variance in the corresponding expression. Additionally to the 
characterization of the most differentially expressed mRNAs and 
their potential functions in relation to SNF, we implemented 
a functional analysis via MapMan for the transcriptomes of 
both tissues. The functional classification of all protein-coding 
chickpea genome sequences will pave the way for future func- 
tional analyses of chickpea mRNAs, also retrospectively (as in the 
present study) and not only from root or nodule, but from any 
tissue of interest. 
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